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Abstract— Canonical Polyadic (or CANDECOMP/PARAFAC, 
CP) decompositions (CPD) are widely applied to analyze high 
order tensors. Existing CPD methods use alternating least square 
(ALS) iterations and hence need to unfold tensors to each of the 
N modes frequently, which is one major bottleneck of efficiency 
for large-scale data and especially when N is large. To overcome 
this problem, in this paper we proposed a new CPD method 
which converts the original iVth (N > 3) order tensor to a 3rd- 
order tensor first. Then the full CPD is realized by decomposing 
this mode reduced tensor followed by a Khatri-Rao product 
projection procedure. This way is quite efficient as unfolding to 
each of the N modes are avoided, and dimensionality reduction 
can also be easily incorporated to further improve the efficiency. 
We show that, under mild conditions, any iVth-order CPD can 
be converted into a 3rd-order case but without destroying the 
essential uniqueness, and theoretically gives the same results as 
direct iV-way CPD methods. Simulations show that, compared 
with state-of-the-art CPD methods, the proposed method is more 
efficient and escape from local solutions more easily. 

Index Terms — CP (PARAFAC) decompositions, tensor decom- 
positions, Khatri-Rao product, alternating least square. 



I. Introduction 

Higher-order tensors (multi-way arrays) have gained in- 
creasing importance as they are often more natural representa- 
tions of multi-dimensional data than matrices in many practical 
applications. As one of the most fundamental problem in 
tensor data analysis, tensor decompositions attempt to finding 
informative representations (e.g. dense/sparse, low-rank repre- 
sentation) of multi-dimensional tensors data. Tensor decom- 
positions are very attractive and versatile because they take 
into account such as spatial, temporal and spectral information, 
and provide links among the various extracted factors or latent 
variables with desired physical or physiological meaning and 
interpretation (TJ, (2). 

As one of the most important tensor decomposition mod- 
els, Canonical Polyadic (CP), also named as CANDE- 
COMP/PARAFAC decomposition (5J, (4| has been extensively 
studied and found many practical applications. One big advan- 
tage of CPD is that the factors are essentially unique under 
mild conditions, which makes it very useful in the cases 
even when no or only very limited a priori knowledge is 
available on the factors. In the CP model, the matricizations 
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(unfoldings) of a tensor are just the products of one factor 
matrix with another matrix which is the Khatri-Rao product 
of all the remaining factors. This has led to the widely adopted 
alternating least square (ALS) methods for CP decompositions, 
e.g. see (5). Unfortunately, in this way we have to unfold the 
tensor to each of the total N modes frequently, which is one 
major bottleneck of efficiency of CPD methods. 

Our recent results showed that, for CP decompositions, 
once only one factor has been correctly estimated, all the 
other factors can be computed uniquely and efficiently by 
using a series of singular value decomposition (SVD) of rank- 
1 matrice^] no matter whether there exist some collinear 
components in factors (6). This motivated us to perform blind 
source separation (BSS) on one single mode first and then 
use efficient rank-1 approximation methods to recover the 
other factors, which has led to the CP-SMBSS method for 
CP decompositions. The CP-SMBSS is quite useful if some 
a priori knowledge about components in at least one mode is 
available. In [ 7 ] independence of components was incorporated 
into CPD. In this paper we deal with the case where such a 
priori information is completely unavailable. 

The following notations will be adopted. Bold capitals (e.g., 
A) and bold lowercase letters (e.g., y) denote matrices and 
vectors, respectively. Calligraphic bold capitals, e.g. ^, denote 
tensors. Unfolding (matricization, flattening) of a tensor ^ G 
r /ix/ 2 x-x/jv i n m ode-n is denoted as Y (n) G R / - xn p#- / p, 
which consists of arranging all possible mode-n tubes (vectors) 
as the columns of a matrix [2]. denotes the Khatri-Rao 
product (column-wise Kronecker product) of matrices and 
OfcL fcl AW = aWqA^-OA^OAW with 
fci > &2- Readers are referred to (TJ, (2) for the notations and 
tensor operations. 

II. CP Decompositions Based on Model Reduction 

A. Review of CP Decompositions 

A CP decomposition (factorization) of a tensor ^G 
^hxi 2 -xi N can be f ormu i a ted as 

y = £A J a«oaf>...oaf>+£, (1) 

J = l 

where matrices A( n ) = [a^ n) , a^ n) , • • • , a^ n) ] G R InXj , 
n G M = {1,2, ••• , N} 9 consists of unknown components 
8Lj n \ e.g., latent source signals, o denotes the outer product, 
and £ is the fitting error, see Fig [I] for the illustration of 

l A matrix Y is rank-1 if and only if Y = uv T , where u and v are any 
nonzero vectors. 
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Fig. 1 : Illustration of CP decompositions of a 3rd-order tensor 
y G ]R / i x/ 2 x/ 3 (ignored the noise), where the factors A^ n ^ = 
(n) a^ n) • • • a^ n) ] G 1 
as their columns, n 
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/nXj contain the latent component 

= 1,2,3. 



CPD of a 3rd-order tensor. From ([I]), the multiway tensor 
is represented as a linear combination of outer products of 
vectors (i.e., rank one tensors), which can be regarded as a 
generalization of matrix singular value decomposition (SVD) 
to tensors [2]. As the scalar factors Aj can be absorbed into a 
factor matrix A^) by letting = \ja^\ Vj, we also use 
y « [A^, A( 2 \ • • • , AW] as a shorthand notation of (l). 

To solve CPD problems, alternating least square (ALS) 
methods are widely employed. Consider the mode-n matri- 
cization of ^: 



(n) 



A (n) B (n)T ? ( n= ! 5 2,... , N.) 



where 



B^ = O A 



(k) 



(2) 



(3) 



In standard ALS methods, factor matrices A^ n ) are updated 
as Y( n )[B( n ) T ]t alternatively for n = 1,2, ••• , TV, where t 
denotes the Moore-Penrose pseudo inverse of a matrix. As the 
matrix B^ n ) is often quite huge, some tricks were proposed 
to simplify the computation of Y^[B^ T }^ (see |2|). 

B. CP Decomposition Based on Mode Reduction 

From the above analysis, when the number of modes N is 
large, the ALS methods often suffer very slow convergence 
speed as they have to unfold the tensor to each mode matrix 
frequently. A number of authors have made efforts to improve 
the efficiency of CPD algorithms, e.g. see (8), (5). In this 
paper we consider a new way to conquer this problem, i.e., 
reducing the number of modes to accelerate the convergence of 
CPD algorithms. As preliminary we need the following tensor 
unfolding operation. 

Tensor transposition (To). Given a CP tensor ^ as ([T]), the 
transposition of ^ is a tensor of I Pl x I P2 • • • x I PN obtained 
by exchanging the roles of A^ n ^ in G AT) accordingly. In 
other words, tensor transposition re-permutes the order of 
dimensions of y. For example, [A^ 1 ), A^ 2 \ . . . , A^)] is 
a tanspose of ^, where (ii, ^2, . . . , ijv) is a permutation of 
(1, 2, . . . , TV). Different to matrix cases, there are an exponen- 
tial number of ways to transpose an TVth-order tensor (lOj. 

Tensor Unfolding. Consider the vectorized version of ([T]) 
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the Khatro-Rao product, in (HJ we can represent some Khatri- 
Rao products of sequential a^ by new vectors , i.e., 



bf) 0b f- 
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and the decreasing sequence n^ (k = 1, 2, . . . , K) denotes a 
split of (TV, TV - 1, . . . , 1) with n K = N and 7V = 0. Rewrite 
([5} in the form of ([I]) and we obtain the reshaped tensor 

j 
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y and y^ K ^ have exactly the same entries but they are 
arranged in different order. Particularly, from ^ there holds 
that 



B w = Q nfc 



i+i 



p=n k 



(8) 



In other words, the tensor unfolding operation actually groups 
and replaces the original factors by their Khatri-Rao products. 
For example, [AW©^ 2 A«, A^" 1 ) ©AW] is a 3-way 
unfolding of ^ with the dimensionality of I\ x (n^T 2 2 Ik) x 
(In-Jn); and [A^, . . . , A^" 2 ), A^" 1 ) AW] is an 
(N — l)-way unfolding. For simplicity we also use a notation 

y{N-li = y{l,2,...,N-2,(N-l)QN} tQ ^ ^ unfold . 

ing denoting that the last two modes are combined. 

By using the tensor transposition and unfolding, the factors 
j±W can be arbitrarily grouped by their Khatri-Rao products, 
thereby leading to mode reduced new tenors. 

It is known from [6], once B^ (k G JC) have been 
estimated, the original factors A^ n ^ in G AT) can be estimated 
immediately and essentially uniquely by using the Khatri-Rao 
product projection procedure, thanks to the special structure of 
B( fc ). Khatri-Rao projection generally can be computed very 
efficiently via, e.g., SVD |6|. A more comprehensive discuss 
on this topic will be detailed in the next section. For simplicity, 
we denote this procedure by KRProj(B^). As 3rd-order 
tensors are the simplest model which has uniqueness guarantee 
under mild conditions, we are interested in converting an TVth- 
order tensor (N > 3) into a 3rd-order tensor in this paper. 
Finally, by using the above tensor transposition/unfolding 
and Khatro-Rao projection, we propose the following mode 
reduction based CPD method (MRCPD): 

Algorithm 1 The General MRCPD Algorithm 

Require: J, and a CP algorithm \I>. 

Let y^^JGW, G (2) , G (3) ] be a 3-way unfolding of y. 
Let(G( 1 ),G( 2 ),G( 3 )) = ^Cy {3} ). 

A( n ) (n G AT) are estimated via efficient Khatri-Rao 
projection KRProj(G^), k = 1, 2, 3. 
return A^ n \n = 1, 2, • • • , TV. 



where y and e are respectively the vectorization of y and £ in 
the proper order of dimensions. Thanks to the associativity of 



Compared with direct iV-way CP decompositions, the MR- 
CPD method only requires 3 factors to be estimated at first, 
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and hence frequently unfolding to each of the total N modes is 
avoided. Moreover, TV-way CPD algorithms may more easily 
stuck into local solutions due to the complexity of the model. 
Finally, some well-designed CPD algorithms only work for 
3rd-order tensors, such as the SWATLD method [11] , etc. 
The MRCPD method makes it possible to apply these methods 
to the CPD of tensors higher than 3 directly. In |T2| we 
have considered the special case where = A~ and 

G ( 2 )=A( 2 ). 

C. Efficient MRCPD Incorporating Dimensionality Reduction 

In (6) we stated that once only one factor has been correctly 
estimated, all the remainders can be essentially uniquely re- 
covered. This motivates us to consider applying 3-way CPD on 

~ {3} 

a reduced tensor y such that one factor, say remains 
unchanged while the size of G^ 2 ^ and G^ 3 ^ is significantly 
reduced: 

. Step 1: Let ^^[G^jG^jG®], where G^ and 
G( 3 ) are factor matrices obtained by reducing the number 
of rows of G^ 2 ) and G^ 3 \ respectively. 

. Step 2: Run 3-way CPD on ^ {3} to obtain G^; 

. Step 3: G( 3 )0G( 2 ) = Y$ T G^ T \ 

• Step 4: Estimate A (n) (n e J\f) by using the Khatri-Rao 
projection recursively on G^ and G^ 3 ^ 0G^. 

The essential uniqueness of 3-way CPD and the full column 
rank of G^ play critical roles in the success of the above 
method. In [6] we considered G^ = where the factor 
A^ 1 ) with full column rank is uniquely estimated by employ- 
ing BSS methods incorporating proper a priori knowledge. 

The above way is quite efficient because we can signifi- 
cantly reduce the size in two modes. Unfortunately it is not 
always the case that G^ is of full column rank. In such a 
case we can reduce the size of only one factor matrix, say 
G^ 3 ), which is generally with the largest size among the three 
factors, while remaining the two others unchanged. After G^ 
and G^ 2 ) have been correctly estimated, we have 

G( 3 )=Y ( { J(G( 2 )0GW) Tt (9) 

In the next step, all the factor matrices A^ n ) (n £ AO can be 
recovered from G^ k \ k = 1,2,3. Note that the full column 
rank of G^ 2 ^ 0G^ is a necessary condition of uniqueness 
for 3-way CPD (13). Consequently, if the 3-way CPD is 
essentially unique, G^ 2 ^ 0G^ is always of full column rank. 
And hence, we can always reduce the size at lease in one mode 
in practice. 

Now we discuss how to reduce the size of one mode 
of y^ 3 \ Suppose that we want to reduce the size of G^ 3 ^ 
without changing the other factors, we consider its mode- 3 
matricization 

Yg> = G< 3 >(G< 2 >©G< 1 >) T , (10) 

from which we observe that reducing the rows of G^ 3 ^ is 
equivalent to reduce the rows of Y^ J . Hence, the following 
dimensionality reduction techniques may be employed: 



1) Truncated SVD. Let Y ( { 3 3 } } = UDV T be the truncated 

SVD of where D £ R JxJ is a diagonal matrix 

whose diagonal elements consist of the leading J sin- 
gular values. Then Y^ is updated as V by letting 
q(3) =d -i u tq(3) e R JxJ ? w hereby leading to the 
significantly reduced size of 

2) Fiber Sampling. Sometimes we need to maintain the 
physical meaning of the original data, e.g., nonneg- 
ativit>r] In this case we can achieve dimension ality 
reduction by sampling the rows of the matrix Ygf 1151, 
which is just equivalent to sampling the rows of G&>. 

It is worth noticing that the above techniques also provide 
an efficient CPD method for 3rd-order tensors incorporating 
dimensionality reduction techniques. Very often one mode, say 
G^ 3 ) = A^ 3 \ can be of extremely large size. We can reduce 
the size of A^ 3 ^ first and then estimate A^ 3 ^ from This 
way provides a trade-off between accuracy and efficiency and 
is quite efficient for large scale problems. 

D. Issue of Uniqueness 

By model reduction the TV- way CPD is converted into a 
lower-order CPD. The first question is whether the simplified 
model is able to give the consistent results with direct CPD of 
the original TV-way tensor. This is answered by the following 
proposition: 

Proposition 1: Let ^S K ^ be a tensor unfolding of an TVth- 
order tensor 3 < K < N. If both y {K} and ^ have 
essentially unique CPD, the MRCPD is able to give essentially 
the same results as the direct TV-way CPD of 

The proof of Proposition [T] is straightforward from the fact 
that the corresponding Khatri-Rao products of A^ n ^ form a 
solution of y^ K \ from the unfolding procedure. Obviously, 
once the CPD of original TVth-order tensor is not essential 
unique, the corresponding lower-order CPD is certainly not 
unique because each solution of the TV-way CPD can always 
form a solution to the lower-order CPD. From Proposition 
[T] the key point is whether the lower-order CPD is also 
essentially unique if the original iV-way CPD is. If it is not 
true, the MRCPD method may lead to very poor performance 
even if the original TV- way CPD is unique. In the following this 
important issue will be investigated based on the well-known 
uniqueness condition given by Kruskal in 1977 for 3-way 
tensors and then extended for TV-way tensors by Sidiropoulos 
and Bro in 2000 (16), (17): 

Lemma 1 [TtJ : For y=[A (1 U (2 V" ,A (iV) ], if 

EN 
kr A( n) >2J+(7V-1), (11) 
n=l 

then the decomposition is essentially unique, where kr A (n) is 
the Kruskal rank (see (l6j) of A^ n \ n = 1, 2, . . . , TV. 

Obviously, higher-order tensors are more likely to have 
unique CPD than lower-order ones (we always assume that 

2 However, it does not mean that the first method cannot be used for 
nonnegative data analysis. Due to the essentially uniqueness of CPD, the 
resulting factors A^ 1 ) and A^ 2 ) are essentially unique and hence can be 
nonnegative after adjusting the signs of their columns accordingly, no matter 
whether ^ 3 J" and G^ 3 ) are negative or not. See also [14J for related 
discussion. 
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kr(A^) > 2 in CP decompositions as it is a necessary 
condition for uniqueness (13). Hence the left hand side of 
Equation ([TT]) increases no slower than the right hand side 



when N increases.) As a result, even the original iVth-order 
tensor meets the uniqueness condition ( [TT] ), it is still possible 
that the corresponding lower-order CPD is not essentially 
unique. This problem will be discussed below. 

Without loss of generality, for a given tensor 
y=J A(1) » A(2) , A (A °1 we hereafter assume that 
k^Aa) > kr A ( 2 ) > • • • > kr A (Ao > 2 (Otherwise we 
transpose the tensor till this condition is satisfied). 

Proposition 2: Given an TV- way (N > 4) tensor 
y=[A (1) , A (2) , • • • ,A W ], if the uniqueness condition (II) 
is satisfied for then the condition ( [TT] ) is also satis- 
fied for the (N — l)-way unfolding tensor ^ N ~ x ^ = 
[A^, . . . , A^- 2 ), A^- 1 ) AW]. Consequently the CPD 
£ \fl{N-i} - s a j SQ essen tially unique. 
(The proof can be found in Appendix.) 

Corollary: Under the assumptions of Proposition [2] for any 
3 < K < N — 1 there exists at least one K-way unfol ding 
y^ K ^ such that its corresponding uniqueness condition ( |ll) 
holds. 

In this paper we only consider K = 3. It can be 
seen that Proposition [2] also provides a way to unfold a 
given A/th-order tensor to a 3rd-order tensor where the 
uniqueness is guaranteed. For example, if N = 5, J — 
18, and assuming that the corresponding Kruskal ranks of 
the factors are 10,9,. . . ,6, respectively, we can construct 

y{3} = | A (l)^ A (2) 0A (3) 5A (4) 0A (5)j which wm have es _ 

sentially unique CPD. In other words, the Kruskal rank of the 
new factor matrices should be mostly balanced as possible. 
This result shows that the method proposed in [ 12] may fail if 
all the factors are ill-conditioned. In practice, the true Kruskal 
ranks of factor matrices are unknown. We may use the mode 
ranks to decide the optimal way for tensor unfolding. Note 
also that, besides the Kruskal uniqueness condition, in (18) 
some relaxed uniqueness conditions were proposed for 3 -way 
CPD. This will simplify the validation of uniqueness. 

Based on the above analysis, we may find another interest- 
ing feature of the MRCPD method. For a TVth-order tensor, 
in practice we may only be interested in one factor which 
plays a role in the TV-way CP structure and will be used for 
further data analysis tasks such as clustering, classification, 
etc. In this case we may consider a 3rd-order CPD first and 
then the desired factor can be extracted by using the Khatri- 
Rao projection procedure. In this way although we actually 
do not perform a full CPD of the original TVth-order tensor, 
the extracted factor is simply consistent with the one obtained 
from full CPD, according to Proposition [T] 

III. Khatri-Rao Product Projection 

A general optimization problem to perform Khatri-Rao 
projection can be formulated as 



mm 

A( fc ),fce/c 



IH-A^OA^©-.^ 1 )^, (12) 



where JC = {1, 2, . . . , K}. In (12), a given matrix H is 
projected onto a Khatri-Rao product space, and the columns 



of A^ (k G JC) can be estimated independently by solving a 
total of J least square problems 



mm 



,keJC 



aWoaf-^O-.-OaWlH, 



(k) 

where j = 1, 2, . . . , J, hj and aj- ; are the jth columns of H 
and A^ k \ respectively. The solutions of {12} are generally 
unique and can be solved, for example, by the procedure 
described in [6|. On the other hand, from ( fT3l ), we can reshape 
such that H^') « a!f\a!f : ~ 1) • • • a*) T , that is, H«) 
can be considered as the mode-if unfolding of a tensor < K^\ 



(13) 



and ( 13 ) is equivalent to 



mm 



(1) (2) 
3 i 03 i °' 



(k) 

In other words, the optimal aj- (k G JC) can be obtained by 
seeking the rank-1 CPD of the Although this can be 

done with resorting to any standard CPD methods, below we 
consider two relatively simple yet efficient implementations. 

1 ) Parallel Extraction: Consider the mode-n unfolding of 
JC^ and (14) is equivalent to 



•oaf) |||. 



(14) 



mm 



neJC 



(15) 



where v = /c _^ n a^. Hence the optimal a^ is just the 
left singular vector associated with the largest singular value 
of h|^. In this way all the columns of A^ k \ k G JC, are 
estimated uniquely and parallely. Consequently, this way may 
benefit much from parallel computations. 

In order to impose specific constraints on the components 
we often employ the power iterations 



following by a projection operation V respectively 



af^PfaH, v = 7>(v). 



An) 



(16) 



(17) 



For nonnegative constraints V is element-wisely defined as 



V+(x) = max(x, 0), 
or for sparsity constraints 

Vs(x) = sign(x)(\x\ - A), 



(18) 



(19) 



where A is a nonnegative parameter. We repeat (16) and (17) 
till convergence. 

2 ) Tensorial Power Iterations: Analogy to the power itera- 
tion method in the matrix case, tensorial power iterations can 
be derived straightforwardly by using tensor operations. From 
( fT4] ) we have 



(n) 



=^x ia ^ T x 2 a (2)T 



Xn+i a 



'3 

(n+l)T 



(n— 1)T 



x K a. 



(K) T 



(20) 



— x n a^- 



(n)T 



5 



Note that v T v = Yl p ^ n a^ T a^ T . Then the general tensorial 
power iteration is 



(n) X n a j 



aJ n ^P(a) n 0, (21) 

lp/n ll a f y " Z 

where P can be ( [T8] ) or ([19]) to impose desired constraints. We 
repeat ^2A) alternatively for n = 1, 2, • • • , K till convergence. 
This is actually an optimal rank-1 CPD of tensors without 
involving matrix inverse operations and is the extension of the 
power iteration method in tensor scenarios. 

By repeating the above procedure for j = 1, 2, . . . , J all the 
columns of A^ n ) in G JC) can be obtained, which realizes the 
Khatri-Rao projection of H and is denoted by KRProj(H). 

IV. Simulations 

Two performance indices (PI) were used to evaluate the 
performance of the proposed algorithm. The first one is the 
signal-to-interference ratio (SIR), which is defined by 



SIR(a,a) = 101og 10 



(22) 



where a, a are normalized random variables with zero mean 
and unit variance, and a is an estimate of a. The value of 
SIR reflects how well the estimated component (source) match 
the true original one. The second PI measures the fit of the 
estimated tensor to the original tensor which is defined as 



Fit(y,y) = i 



(23) 



where ^ is an estimate of Fitfy,^) = 1 if and only if 
y = y. For synthetic data, ^ is the original noiseless data in 
order to evaluate how robust of the proposed method in respect 
to additive noise. All the experiments were done in MATLAB 
2008a on a computer with Intel 7i 3.33GHz CPU and 24GB 
memory running Windows 7. 

Simulation 1: We generated a 5th-order tensor y using 
the CP model. The elements of each factor matrix A^ n ) G 
R 20x48 were drawn from independent standard normal dis- 
tributions, which means that each factor is undetermined and 
5] n rank(A( n )) = 100 = 2 J + (N - 1). This setting makes 
the problem rather difficult as it is on the boundary of the 
uniqueness condition ( pT) . Finally independent Gaussian noise 
with SNR=20dB was added to the observation tensor. The 
proposed method was compared with the standard CP method 
based on ALS iterations (CP-ALS) in |5], the CP-ALS com- 
bined with line search (CP-ALSLS) (19), and the PARAFAC 
algorithm included in the TV-way tensor toolbox for MATLAB 
|20) (nPARAFAC, ver. 3.20). For these methods the maximum 
iteration number was set to 100. In the MRCPD method, we 
constructed a 3rd-order tensor ^ {3} = t^ 1 * 2 © 3 * 4 © 5 ^ an d 
performed PC A on the mode- 3 of Then we used the 

nPARAFAC method to perform 3-way CPD. Finally all the 
factors were recovered by using tensorial power iterations. 
Their performance over 50 Mente Carlo runs was detailed 
in TABLE [TJ where the Global Convergence Rate (GCR) 
evaluates the ability of escaping from local solutions. The 
corresponding Fit is plotted in Fig|2j which shows that only 



TABLE I: Performance comparison between the algorithms 
in terms of Runtime, Fit, SIR, and GCR averaged over 50 
Monte-Carlo runs. 



Algorithm 


Runtime(s) 


Fit 


SIR 


GCR 


CP-ALS 


59.2 


0.94 


52.5 


34% 


nPARAFAC 


236.5 


0.94 


52.2 


30% 


CP-ALSLS 


59.7 


0.97 


53.1 


66% 


MRCPD 


1.4 


1.00 


54.9 


100% 




0.85 



0.75 



Fig. 2: Fit of each algorithm in 50 Monte-Carlo runs. The 
MRCPD method can escape from local solutions more easily. 

TABLE II: Performance comparison of the four CPD methods 
when they were applied to real image data clustering over 20 
Monte-Carlo runs. 



Algorithm 


Fit 


Runtime (s) 


Accuracy 


CP-ALS 


0.67 


564.2 


63.8% 


HALS 


0.65 


506.8 


65.5% 


nPARAFAC 


0.67 


746.5 


65.4 % 


MRCP 


0.64 


44.6 


67.0 % 



MRCPD found nearly the optimal Fit in all runs. From the 
simulation results, the MRCPD method is much more efficient 
and escape from local minima more easily than the other 
compared methods. 

Simulation 2: We applied the proposed methods to real 
image data analysis, namely the COIL- 100 database pl| . The 
COIL- 100 database consists of 7200 color images of 100 
objects, and 72 images per object which were taken from 
72 different angles. For simplicity, we selected the first 20 
objects for test, and each image was scaled with the size of 
128 x 128. Then a tensor ^ with the size of 128 x 128 x 3 x 1440 
was generated. We set rank J = 10 for all methods. In the 
MRCPD, y {3} = yi 1 ? 4 ? 2 © 3 } and we randomly sampled 100 
fibers in mode-3 in each run. We used the HALS method (22) 
to perform 3-way CPD. For the other methods the maximum 
iteration number was set to 100. Finally, the factor A^ 4 ^ was 
used as features to cluster the original images. As iiT-means 
is prone to be influenced by initial centers of clusters, we 
replicated if -means 20 times for each method. See TABLE 
|Il|for their performance over 20 Monte Carlo runs. From the 
table, we see that the MRCPD method achieved comparable 
clustering accuracy but is significantly faster than the other 
methods. 
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V. Conclusion 
Existing CPD methods are often based on alternating least 
square (ALS) iterations. In ALS methods, we have to unfold 
the observation tensor frequently, which are the major bottle- 
necks of high efficiency. To overcome this problem, in this 
paper we proposed the concept of mode reduction for the 
first time and based on it we developed a new method to 
perform CPD of higher order tensors with TV > 3, efficiently 
incorporating dimensionality reduction techniques. In this way, 
frequently unfolding with respect to each mode is avoided 
and the new method can escape from local solutions more 
easily due to the significantly reduced complexity of model. 
Moreover, a full TV-way CPD may be simply unnecessary if 
only partial factors are desired, but without loss of structural 
information of data. The essential uniqueness of the mode 
reduced tensor was also theoretically investigated. Simulations 
confirmed the efficiency and validity of the proposed method. 

Appendix A 
Proof of Proposition 1 

Before proof we introduce the following Lemma: 
Lemma 2 (Lemma 3.3 in (T3j): Consider matrices A G 

R /ixJ and B e R / 2 xJ if kr A > 1 and kr B > 1, then 

krAQB > min(kr A + kr B - 1, J). 



We have assume that kr A (i) > kr A 



(2) 



> • • 



> kr 



AW 



> 2 



as it is a necessary condition for uniqueness (T3J. Conse- 
quently kr A (n 1 ) A (n 2 ) > min(kr A (n 1 ) + kr A (n 2 ) — 1, J) holds 
for any ni =^ ri2, from Lemma 2. 

Proof of Proposition |2| Consider that y^ N ~^ = 
[A^, . . . , A^- 2 ), A^- 1 ) AW] where the two factors 
with the minimum Kruskal ranks are merged into one mode. 
Let t = kr A (jv-i) A (jv). We only need to show that 
J2k=i ! ^ r A( fc ) + t > 2 J + (N — 2), a sufficient uniqueness 
condition for the CPD of \) {N ~ 1} . 

From Lemma 2, t > min(kr A (Ar-i) + kr A (jv) — 1, J). There 
are two possibilities: 

1) kr A ( at — i ) + kr A( iv) — 1 > J, i.e. t = J and 

2kr A (7v-i) > kr A (7v-i) + kr A (iv) > J + 2. (24) 



N-2 



Note that £ n= i 
t = J, we have 



kr A( n) > (N - 2)kr A( iv-i). From (24) and 



.AT-2 



>(N - 
>\(N 



2) 



J 



that is, Y,k=i krAffc) 



2 

4) J > 0, 
t > 2 J 



J 



[2J- 
- [2J- 



(N- 
- (N - 



2)] 
-2)] 



(25) 



2) kr A (jv-i) + kr A (Ar) 
kr A (7v) — 1. Moreover, 



(N-2). 
1 < J, and hence t > kr A (jv-i) 



> 



kr A (n) 
kr A( n) - 



(kr A(J v-i) +kr A 



(N) 



1) 



(26) 



= E 
>2J + (A/--2). 
The proof is complete. 
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